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We obtained the exact solution of a probabilistic cellular automaton related to 
the diagonal-to-diagonal transfer matrix of the six-vertex model on a square lattice. 

The model describes the flow of ants (or particles), traveling on an one-dimensional 
lattice whose sites are small craters containing sleeping or awake ants (two kinds of 
particles). We found the Bethe ansatz equations and the spectral gap for the time- 
evolution operator of the cellular automaton. From the spectral gap we show that in 
the asymmetric case it belongs to the Kardar-Parisi-Zhang (KPZ) universality class, 
exhibiting a dynamical critical exponent value z = This result is also obtained 
from a direct Monte Carlo simulation, by evaluating the lattice-size dependence of 
the decay time to the stationary state. 

I. INTRODUCTION 

The six-vertex model was introduced in 1931 by Pauling in order to explain the residual 
entropy of the ice at zero temperature. The model turns out to be of great interest for 
physicists and mathematicians of many-body systems due to its exact integrability [Tj. The 
row-to-row transfer matrix of the six-vertex model is the generating function for an infinite 
set of commuting non-trivial charges in involution j2i|. The anisotropic Heisenberg chain, 
or the so called XXZ quantum chain, is one of these conserved charges. Actually, a quan¬ 
tum system is integrable whenever its Hamiltonian belongs to an infinite set of commuting 
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operators. The exact integrability of the XXZ quantum chain is then a consequence of the 
infinite number of commuting charges generated by the six-vertex model. For this reason the 
six-vertex model is considered as a paradigm of exact integrability in statistical mechanics 

On the other hand, the representation of interacting stochastic particle dynamics in terms 
of quantum spin systems produced interesting and fruitful interchanges among the area of 
equilibrium and non-equilibrium statistical mechanics. The connection among these areas 
follows from the similarity between the master equation describing the time-fluctuations on 
non-equilibrium stochastic problems and the quantum fluctuations of equilibrium quantum 
spin chains [THUS EDH2S]. The simplest example is the asymmetric diffusion of hard-core 
particles on the one dimensional lattice (see pjEJES, 28J for reviews), where the time fluctu¬ 
ations are governed by the time evolution operator that coincides with the exact integrable 
XXZ quantum chain in its ferromagnetically ordered regime. 

An important consequence of the above mentioned mathematical connection between 
quantum chains and interacting stochastic problems is that, unlike the area of non¬ 
equilibrium interacting systems where very few models are fully solvable, there exists a 
huge family of quantum chains appearing in equilibrium problems that are integrable by the 
Bethe ansatz on its several formulations (see [5l f29ff3T] for reviews). In these cases, the Bethe 
ansatz enable us to obtain exactly the complete spectrum of the master equation describing 
the time-fluctuations on the non-equilibrium stochastic problem. From the spectrum we can 
compute important properties of the system. As for example, the relaxation time to the 
stationary state which depends on the system size L. It satisfies a scaling relation ry ~ L z , 
where the dynamical exponent z can be obtained from the finite-size dependence of the real 
part of spectral gap between of the two leading eigenvalues ruling the long time regime. 

In addition, although less explored in the literature, there is also a non-trivial connection 
between the transfer matrix operator of two-dimensional vertex models in thermal equilib¬ 
rium and the time evolution operator in discrete time Markov chains [28]. Examples of this 
connection are given by the six-vertex model. It describes a non-local parallel diffusion of 
hard-core particles USED], and also the time-evolution operator of a two steps local cellular 
automaton describing a half-parallel dynamics of asymmetric diffusion of hard-core particles 
[28, 323- Another application is given by the ten-vertex model introduced in [33] to describe 
a parallel dynamics of a traffic ffow stochastic cellular automaton. 
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The connection of the KPZ dynamics with the row-to-row transfer matrix of the six- 
vertex model was first observed by Gwa and Spohn [IB]. This result for the six-vertex model 
on a cylinder was confirmed by several works including the analytical computation from 
Bukman and Shore H3. Furthermore, more recently a rigorous proof of the KPZ exponent 
was also obtained for the six-vertex model on a quadrant with parameters on the stochastic 
line by Borodin, Corwin and Gorin [183 • It is interesting also to mention that in [T9] we can 
find a recent review of theoretical and experimental realizations of the KPZ dynamics. 

In the present work we study a simple one-dimensional probabilistic cellular automaton 
related to the diagonal-to-diagonal transfer matrix of the six-vertex model. This model is 
inspired in the idea of Alcaraz and Bariev [33] and describes the flow of particles, ants, 
automobiles, or some other conserved quantity in a one-dimensional discrete chain. Despite 
the diagonal-to-diagonal six-vertex model being known to be solvable for quite a long time 
[34] . to our knowledge the spectral gap for this model with Boltzmann weights associated 
to stochastic processes was not considered previously. In this work we calculate the spectral 
gap for this model and found that the associated probabilistic cellular automaton belongs 
to the Kardar-Parizi-Zhang (KPZ) universality class with a dynamical exponent z = f [35]. 

The layout of this papers is as follows. In Section II we introduce the probabilist cellular 
automaton and show its correspondence with the row-to-row transfer matrix of the six vertex 
model. In section III we solve exactly the model using a matrix-product ansatz. In section 
IV we calculate numerically the spectral gap related to the highest eigenvalues of the time- 
evolution operator and obtain the dynamical critical exponent. In Section V we perform 
Monte Carlo simulations of the cellular automaton. Finally in section VI we conclude our 
paper with a summary of the main results. 


II. THE STOCHASTIC CELLULAR AUTOMATON AND ITS CONNECTION 
WITH THE EXACTLY INTEGRABLE SIX-VERTEX MODE 

The main ingredients of a cellular automaton are the associated configuration space and 
the dynamical rules defining its time evolution. The cellular automaton we consider is 
defined on a L-site chain where at each site we attach a variable that have four possible 
values, and the configuration space has then 4 L components. We can associate to the four 
possibilities of site occupation arbitrary objects. We chose here a simple formulation of the 
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automation where the four possibilities corresponds to a vacant site, a site occupied by a 
single awake or sleeping ant, and a site double occupied by an awake and a sleepy ant. We 
should stress here that the probabilist cellular automaton we introduce is not intended to 
describe the realistic motion of real ants. For a more realistic motion we should also consider 
the effect of the pheromones in the ants motion, like e.g., in [36]. 

We consider n ants traveling in a one dimensional tunnel that links L small craters (sites), 
as shown in Fig. [l] On the craters the ants are on a sleeping or awake state. The craters 
are so small that at most it supports two ants and no more than one awake ant. The four 
possible occupations of a given crater are then: no ants, one ant (sleeping or awake) or two 
ants (one sleeping and another awake). 

The dynamical rules are defined as follows. At a given unit of time all the ants travel to 
the next nearest crater to the left, or stay in the same crater. The travel happens with the 
following rules: 

a) If the ant is awake and alone in the crater it sleeps with probability p or travels to the 
left crater with the probability 1 — p. 

b) If when the ant arrives to a crater it finds another ant sleeping, since they do not want 
a dispute (the connecting channels among the craters only allow a single traveling ant) 
with probability 1 the ant travels to the left leaving the other one sleeping. 

c) If in a given crater there is a sleeping ant with a probability q it awakes and travels 
to the left and with probability 1 — q stay sleeping. 



FIG. 1: An illustrative example of a one dimensional set of tunnels connecting small craters. 
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The allowed motions, in a unit of time, can be represented by the arrow configurations of 
the six vertices shown in Fig. [2j The motions in Fig. [2] are from the top to the bottom of 
the figure. A vertical down arrow (J,) represents an ant sleeping while the diagonal arrow 
{y/') an awake ant, and lines with no arrows indicate the absence of ants. For simplicity, 
we consider periodic boundary condition. The conservation of the number of ants is then 
a consequence of the ice type rules defining the vertices, that demand an equal number of 
outgoing and incoming arrows in a given site. 


a o 


a x 


b i 


b 2 c x c 2 








FIG. 2: The fugacities of the six allowed vertex configurations of the equivalent vertex model. 
They correspond to the probabilities of motion in the probabilistic cellular automaton. 

Let us denote as ay the possible occupations of the crater j by oij = 0,1, 2 and 3, that 
indicates that there is no ant, a sleeping ant, an awake ant, and a sleeping and awake 
and, respectively. The probability distribution P({a},f) of finding the system in the state 
|{a}) = |ai, a 2 ,..., a v ..., a^) at time t is given by the components of the probability vector. 

|J’({«},«)) = E J ’({“}- t )l{“}>' (!) 

M 

where the summation is performed over all possible configurations {a} = {ay, ...,cul} satis¬ 
fying the ice type rules. The time evolution of the stochastic cellular automaton is given by 
the master equation: 

\P({a},t + At) = T D . D \P({a},t), (2) 

where the transition matrix T D _ D corresponds to the diagonal-to-diagonal transfer matrix of 
a six-vertex model with fugacities a 0 , a x , b x , b 2 , c x , c 2 related, respectively, to the probabilities 
1,1,1 —p, 1 — q, p, q of the cellular automaton (see Fig. [2]). The components ({a}\T D _ d\W}) 
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give the probability of motion from configuration {a} to {a'} and are given by the product 
of the vertex fugacities connecting the configurations: 

({a}\T D . D \{a'}) = (3) 

where n 0 , ni, n 2 , n 3 ,n 4 , n 5 are the numbers of vertices with fugacities a 0 , cq, b\, b 2 , Ci, c 2 , re¬ 
spectively (see Fig. [3]). 

oij a 2 a 3 a 4 a 5 a z-i 

{a}), t 

{«'}), t+\t 

a\ a ' 2 a ' 3 a ' 4 a ' 5 a' L _ l a' L a\ 



FIG. 3: Pictoric representation of the diagonal-to-diagonal transfer matrix. The variables ay and 
a'j specify the possible arrow configurations at site position j at times t and t + At, respectively. 

The spectral properties of the six-vertex model are characterized by the parameter [5]: 

_ &1&2 — C1C2 + aiao , 

~ 2 \^b 2 ’ 1 ’ 

that for the cellular automaton is given by 

1 _ (p+ q) 

A = (5) 

1 -q 

III. DIAGONALIZATION OF THE TRANSFER MATRIX 

In this section we solve the eigenvalue equation for the diagonal-to-diagonal transfer 
matrix Td~d introduced in the previous section. As a consequence of the conservation of 
arrows, and due to the periodic boundary condition, the transfer matrix Td~d can be split 
into block disjoint sectors labeled by the number n (n = 0,1,..., 2 L) of arrows and by 
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the momentum P = ^-j (j = 0,1 ,... ,L). We want to solve, in each of these sectors, the 
eigenvalue equation 


^n,P | ^n,P >~ T D _ D \Vn,P >, 


( 6 ) 


where A nP and |\I/ njP > are the eigenvalues and eigenvectors of T D _ D , respectively. These 
eigenvectors can be written in general as 

(*) 

\Vn,p >= Y] Yl a n foi, ■ ■ ■ > x n ) [xi, Q!i;...; x n , a n >, (7) 

{*} {«} 

where <j>^ (aq, ■ ■ •, x n) is the amplitude corresponding to the arrows configuration where 
n arrows of type aq,..., a n are located at positions x\, ..., x n , respectively (where a 3 = 1 
correspond to a vertical arrow and a 3 = 2 to a diagonal arrow). Finally, the symbol (*) in 
the last equation means that the sums in {x} and {a} are restricted to the sets obeying the 
ice rule. 

Since |\1 / n ,p > is also an eigenvector of the translation operator with momentum P = yj 
(j = 0,..., L — 1), the amplitudes satisfy 

■ ■ ■ i Xn ) _ -iP /o', 

<,...^(^l + l,...^n + l) " ’ U 

for x n < L — 1, whilst for x n = L 


_ ^ai,...,Qn( Zl ’ • • • ’ L) _ _ iP /g\ 

The solution of the eigenvalue equation ([6]) can be obtained by an appropriate ansatz 
for the unknown amplitudes 0 q, ( x >, ■ • • i x n)- Although the model can be solved by the 

coordinate Bethe ansatz EU. we are going to formulate a new matrix product ansatz (MPA) 
[B8l [40] due its simplicity and unifying implementation for arbitrary systems. This new MPA 
introduced in J38l 14(1] can be seen as a matrix product formulation of the coordinate Bethe 
ansatz and it is suited to describe all eigenvectors of integrable models, including spin 
chains pTSl fTTTj . stochastic models jTTl 1TT| and transfer matrices [Hl - HBj . According to this 
ansatz, there is a correspondence between the amplitudes of the eigenvectors and matrix 
products among matrices obeying special algebraic relations. In the present work, in order 
to formulate the MPA we make a one-to-one correspondence between the configurations 
of arrows and the products of matrices. The matrix product associated to a given arrow 




configuration is obtained by associating a matrix E to the sites with no arrow, a matrix AA') 
to the sites with a single arrow of type a (a = 1,2), and finally the matrix E ~ l for 
sites with two arrows. The unknown amplitudes in (J7|) are obtained by associating them to 
the matrix product ansatz 


b p 

«i, 


1; ... )Qn (aM, ...,X n ) ^i-l J 4(«i) jB x 2 -xi-l A (a 2 ) . . . E x n -x n -i-l A ( an ) E L-x n 

if there are no sites with two arrows, and by associating 

■ ■ ■ ,x n ) ^ e^A^E**-*'- 1 ^ ■ ■ ■ 


( 10 ) 


( 11 ) 


E x i~ x i-^~ l E ~ l E x j+i- x i~ l . . . E ^n-x n -i-l A {a n ) E L-x n 

if there are two arrows at position Xj. The other cases follows straightforwardly. 

Actually E and A^ a > are in general abstract operators with an associative product. A well 
defined eigenvector is obtained, apart from a normalization factor, if all the amplitudes are 
related uniquely, due to the algebraic relations (to be fixed) among the matrices E and A^ a \ 


Equivalently, the correspondences (10) and (11) implies that, in the subset of words (prod¬ 


ucts of matrices A A) and E) there exists only a single independent word (’’normalization 
constant”). The relation between any two words is a c number that gives the ratio between 
the corresponding amplitudes. Moreover, since the eigenvectors have a well defined momen¬ 
tum P = ^l (/ = 0,..., L — 1), the relations (pi) and (^| imply the following constraints for 


( 12 ) 


the matrix products appearing in the ansatz (10) 


E X!-1 A ( ai ) E X2-X!-1 A (ct 2 ) . . . E x n —x n -\-l A {a n ) E L-x n _ 

e ~iP E x i A (ai) E x 2 -xi-l _ _ _ A (a n ) E L-x n -l 

for x n < L — 1, and for x n = L 

E x !-i A (a ± ) E x 2 - Xl -i A {a 2 ) ■ ■ ■ E Xn ~ Xn ~ 1 ~ 1 A^ an ^E L ~ Xn = 

e ~iP A (a n ) E xi-1 A (ai) , . . A (a n - 1 ) E L-Xn- l-l 


(13) 


The relations (12) can be easily solved by identifying the matrices A^ as being composed 


by n spectral dependent matrices A^ j = (1 ,,n), 

n 

A a = Y,K^E, a = 1,2, 

3 =1 

satisfying the algebraic relation 


(14) 


EA, = e ikj Ai,.E. 


( 15 ) 
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By inserting (14) and using (15) into (12) we verify that the spectral parameters kj (j = 


1, n ) are related to the momentum P = (l — 0,..., L — 1) of the eigenvector: 

n 

P = E k r 

3 = 1 


(16) 


On the other hand, by inserting (14) and using (15) into the boundary equations (12) we 
obtain the algebraic relations 


E x 'A kl E x *- x 'A k2 ---E 


L ~xn-i Ak = A k E xl ~ l A k ---A k v E L ~ Xn -\ 

r^n •*'1 'Vn— 1 


(17) 


The eigenvalue equation (|6]) give us two kinds of relations for the amplitudes (10) and 0- 
The first one is related to those amplitudes without multiple occupancy, and the second one 
is related to those amplitudes with multiple occupancy. The first kind of relations, after some 
algebraic manipulations following [44] (with t — 0 in [44], where t represent an interaction 
parameter of the interacting five vertex model), give us the eigenvalues A n ,p(ki ,..., k n ) of 
the transfer matrix 

n 

An,p(&i, • • • ,k n ) = JjAi(%), (18) 

3 =1 

where 


A i(k) = A ^\k) = (b 2 + b x e lk ± \J(b 2 + he^) 2 - - c 2 ci) j 

= \ (l _ Q + (1 ~ ^ e%k ^ \/[l - q + (1 - p)e ik } 2 - 4e ifc (l - p - q)^ 


(19) 


On the other hand, the relations coming from the configurations where two arrows are 
at same position fix the commutations relations among the matrices A k .: 


A kj A kl s(kj, ki)A kl A kj , 


( 20 ) 


where 

/, , \ _ _A\(ki)Ai(kj)bi — Ai(kj)(b 2 bi — c 2 ci + ai) + a\b 2 
1)1 Ai(ki)Ai(kj)bi — Ai(ki)(b 2 bi —<c 2 ci + ai) + a\b 2 

_ Ai(fc;)Ai(fcj)(l - p) - A\{kj)(2 -p-q)+ l- g 
Ai(ki)Ai(kj)(l -p) - Ai(k t )(2 -p - q) + l-q' 


( 21 ) 


The relations with more than two arrows at same positions are automatically satisfied, due 
to the associativity of the algebra of the matrices A k . (j — 1,, n). 
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Finally, the up to now free spectral parameters {ki} are fixed by the nonlinear Bethe 
equations of the diagonal-to-diagonal six-vertex model 


e ikjL _ 


Y[s(kj,ki), j = l,...,n, (22) 

1=1 

that is obtained by imposing the consistence between boundary relations ([9]) and the com¬ 


mutation relations (15) and (21) 


IV. THE SPECTRAL GAP 


In order to complete the solution of any integrable model we need to find the roots of 


the associated spectral parameter equations (Eq. (22) in our case). The solution of those 
equations is in general a quite difficult problem for finite L. However, numerical analysis on 
small lattices allows us to conjecture, for each problem, the particular distributions of roots 
that correspond to the most important eigenvalues in the bulk limit L —» oo. Those are the 
eigenvalues with larger real part in the case of transfer matrix calculations. The equations 
we obtained in the last section, up to our knowledge, were never analyzed previously for 
either finite or infinite values of L in the asymmetric case bi ^ b 2 and C\ c 2 (p q). 
Actually, while a detailed numerical analysis of the Bethe ansatz equations for the row-to- 
row transfer matrix of the six-vertex model, or the XXZ chain, was previously done [471 48] , 
the diagonal-to-diagonal model was previously analyzed only in the symmetric case and with 
A > 1 [g]. 

The corresponding Bethe equations for the diagonal-to-diagonal transfer matrix of the 
six-vertex model are quite different from those of the row-to-row transfer matrix. In order 
to simplify our analysis we are going to restrict ourselves hereafter to the case where p = 1 


(ci = 1 and b\ — 1 — p — 0) in a half-filled chain. In this case we have A = | and, from (19) 
and (21), 

r iki _ + Q~ 1] 

^ J 

q 

and 


S{k j ,k l ) = 


[Ai {kj) - 1] 


l A i( k i) - !]' 

By inserting (23) and (24) in (21), the Bethe equations can be rewritten as 

\L+1 L 


(23) 

(24) 


Ai (kj 


_Ai(%)(Ai(%) + q - 1) 


-1) J 


r 


hia.w -1]= r . 


(25) 


1=1 
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where we set n — L (half-filled chain). If we parametrize Y = —a L e l9L , with a > 0 and 

5) are given by 

1 - y 3 {q - 1) ± y/(yj(q-l) - l) 2 - 4 y 3 


6 e (—f, f), the roots A i(kj) of (25) are given by 


A i(kj) = A = 


y 3 = ae ie e i27TU ~^ /L = 


2 Vi 


Ai (kj) - 1 


(26) 


Ai(%)(Ai(%) + q- 1)’ 


j = A. 


For a given choice (Ai(fcj)} of the above set of roots we have from (25) and (26) 

(- 1) 1 


a L e i9L = 


q‘ 




(27) 


i=i 


We have solved numerically the above equations for several values of q. The stationary 
states A° L = Al,p = 1 eigenvalue is obtained by setting P = kj = 0 and by choosing 


A° = A^ikJA^ih) • • • Ai +) (A; i _ 1 )AS +) (A; L ), 


(28) 


where A i\kj) are given by (26). On the other hand, the eigenvalue having the second 
largest real part belongs to the eigensector with momentum P = ^kj = ^, and is obtained 
by choosing 

Ai = A LtP = AUkJA+ih) ■ ■ ■ A+(fc i _ 1 )Ar(fc L+1 ). (29) 

In Table [I] we display our numerical solutions for |Ai|, with p — 1 and q = 0.7, 0.8, 0.9, 
for several values of L. The eigenvalue with second largest real part (29) determines the 
relaxation time and the dynamical exponent z. By rewritten A^ = |Ai|e* 7 = e~A( L )+*7(0 5 
we have 


fq(L) ~ L~ z . (30) 

In Fig. [d] we show In (f q (L)) versus ln(L) for p — 1 and q = 0.7, 0.8, 0.9. As we can see 
from this figure the energy gap give us a dynamical critical exponent z — |. Consequently, 
the probabilistic cellular automaton introduced in Section II and the related asymmetric 
six-vertex model belongs to the KPZ universality class [35J. 


V. MONTE CARLO SIMULATIONS 

In order to check and illustrate the results of last section we present some Monte Carlo 
simulations (MCS) of the stochastic cellular automaton whose dynamic rules were defined 
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L 

K(L,p=l)\ 

q = 0.7 

q = 0.8 

q = 0.9 

4 

0.94140144138869053 

0.96301968197766430 

0.98238800386894887 

10 

0.98538385114230209 

0.99085409804497404 

0.99568145771034344 

50 

0.99874631037139583 

0.99921749651303016 

0.99962661788291240 

100 

0.99955908669853510 

0.99972484139884665 

0.99986844854358292 

200 

0.99984424258354831 

0.99990280245468144 

0.99995359703973385 

300 

0.99991524776815632 

0.99994711263234282 

0.99997476777514405 


TABLE I: The real part of the second largest eigenvalue of Te>_d for lattice sizes L = 
4,10, 50,100, 200, 300 and with parameter p = 1 and q = 0.7, 0.8,0.9. 



FIG. 4: The logarithm of f q (L), given in (30) as a function of In L for the stochastic model with 
parameters q = 0.7, 0.8, 0.9 in the in the half-filling sector n = L. The dynamical exponent is close 
the expected value z = |. 


in Sec. II. We consider, as in last section, the case where the number of ants is equal to the 
lattice size n = L (half-filled lattice). 

It is not simple to calculate the dynamical critical exponent z of models in the KPZ 
universality class by measuring directly the decay-time of observables starting on a given 
initial condition. A known way is by considering the time-correlations of tagged particles 
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FIG. 5: The ratio a s (t) / a s (oo), as defined in (31), of the standard deviation of the local density of 
sleeping ants, as a function of time. The results were obtained from Monte Carlo simulations for 
the cellular automaton in the sector with n = L (half-filled) ants and with parameters p = 1 and 
q = 0.75. At the initial state all the L craters are occupied by a single sleeping ant. The dashed box 
gives the region used for the linear fit to extract the typical decay time tl (see text).The results 
are the average of 10 7 samples for each lattice size. 


(ants in the present case) |32l EQ] • We verified that, alternatively, this exponent can also be 
calculated by measuring the time evolution of the variance of a local operator in a given 
site. In particular by defining n s (i,t ) as the occupation number (0 or 1) of a sleeping ant at 
site i and time t, we consider the measure 


a‘(t) = 


\ 


|n*(i,«)- < n ° >c°] 2 >, 


(31) 


1=1 


where < n s >oo is the asymptotic average number of sleeping ants in a given site. 

We consider two kinds of initial states. Initially we consider a state not translational 
invariant. We take L /2 consecutive craters double occupied (one sleepy and one awake ant) 
and the remaining craters empty. In Fig. [5] we show the ratio a s (t) / a s (oo) for lattice sizes 
L = 10,20,30,40 and 50, and parameters p — 1, q — 0.75. We estimate from this figure 
the typical times tl where the systems of size L reach the stationary state. These typical 
times are estimated from the crossing with the value 1 of the straight line obtained from 
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FIG. 6: The typical decay time tl, for lattice sizes L = 20,20,30,40 and 50. The dashed box 
gives the region used for the linear fit to extract the typical decay time 77, (see text). The squares 
are obtained from Fig. [5j where p = 1 and q = 0.75. The circles are the results obtained for the 
symmetric case p = q = 1/2. Both figures were obtained by the same procedure and number of 
samples. The estimated values for the dynamical critical exponent z are also close to the values 
3/2 and 2. 

the fitting in the region represented by a dash box in Fig. [5] As we can see the time decay 
tl increases with the lattice size. The dynamical critical exponent z is calculated from the 
finite-size behavior ry ~ L z . In Fig. [6] we fit the values of ln(ry) with ln(L) and obtain the 
estimate z = 1.54 ± 0.05, in agreement with the expected value z = 3/2. In the case where 
the cellular automaton has parameters p — q — 1/2, i. e. the symmetric case, we expect 
that the dynamics is pure diffusive. For the sake of comparison we also show in Fig. [6] the 
lattice dependence of the typical decay time obtained for this symmetric case. We obtain in 
this case z = 2.02 ± 0.05, in agreement with a standard diffusive behavior where z — 2. 

We also consider initial states that are translational invariant, as for example the one 
where we have a single awake ant on each crater (n = L ). We verified that by using the 
stochastic rules that define the cellular automaton introduced in section II it is difficult to 
extract the dynamical critical exponent z. In this case the typical times where the system 
achieve the stationary state is quite small. We can however extract some reliable estimate 



15 



FIG. 7: The ratio a s (t) / a s (oo), as defined in (31), of the standard deviation of the local density 
of sleeping ants, as a function of time. The results were obtained from Monte Carlo simulations 
of the cellular automaton with a small local change of the dynamical rules defined in Sec. II (see 
the text). As in Fig. [ 5 ] the parameters are p = 1 and q = 0.75 and 10 7 samples were considered 
for each lattice size. The dashed box delineate the region where a linear fit is done to extract the 
typical decay time tl (see text). 


of z by introducing a small local change in the dynamical rules, that keeps the model on 
the same universality class of critical behavior. The modified dynamics keeps the number of 
ants conserved but the last crater (L) is special. If we have an awake ant at this crater or 
there is no ants coming to occupy this crater (coming from crater 1) the rules are the ones 
that define the cellular automaton (see Fig. [ 2 ]). However if one ant is coming (from crater 
1) and there is no awake ant at crater L , there is two possibilities with distinct transition 
probabilities. If there is no ant at crater L with probability p the arriving ant stays sleeping 
at the crater L, and with probability 1 — p goes to the crater L — 1, staying awake. If there 
is already a sleeping ant with probability one the awake ant that comes from crater 1, goes 
directly to crater L — 1 staying awake, and leaving at crater L the sleeping ant. In Fig. [7] we 
show the ratio a s (t)/a s (oo ) for the same lattice sizes and parameters as in Fig. [5] In Fig. [8] 
we show the dependence of the typical times tl to reach the stationary state as a function 
of the lattice size L. We obtain from this last curve the estimate z = 1.47 ± 0.05, again 
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FIG. 8: The typical decay time tl, for lattice sizes L = 20, 20,30,40 and 50 obtained from Fig. [7J 
The estimated value for the dynamical critical exponent z is shown in the figure. 

in agreement with the KPZ expected value z = 3/2, as predicted by the exact solution of 
section IV. 


VI. CONCLUSION 

In the present work we obtained the exact solution of a probabilistic cellular automaton 
related to the diagonal-to-diagonal transfer matrix of the asymmetric six-vertex model in 
a square lattice. The model describes the flow of ants (or particles) on a one-dimensional 
lattice where we have sleeping or awake ants on a ring of L craters. The solution was 
obtained by a matrix-product ansatz that can be seen as a matrix product formulation of 
the coordinate Bethe ansatz. Solving numerically the Bethe ansatz equations we calculated 
the spectral gap of the model. From the f ini te-size dependence of the spectral gap we verified 
that the model belongs to the KPZ universality class displaying a dynamical exponent z — |. 
This result was also verified from MCS. The evaluation of this exponent using MCS can be 
done by measuring the time evolution of the standard deviation of local operators. The 
lattice-size dependence of the decay time to the stationary state give us an estimator for z. 
Reliable results are obtained by initiating the system in a non translational invariant state. 
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We can also obtain estimates for z starting with translational invariant states ,provide a 
small local change of the dynamics, playing the rule of a ’’defect”, is done. 
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